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ABSTRACT 



We integrate for the first time the hydrodynamic Hall-Vinen-Bekarevich- 
Khalatnikov equations of motion of a 1 So-paired neutron superfluid in a rotating 
spherical shell, using a pseudospectral collocation algorithm coupled with a time- 
split fractional scheme. Numerical instabilities are smoothed by spectral filtering. 
Three numerical experiments are conducted, with the following results, (i) When 
the inner and outer spheres are put into steady differential rotation, the viscous 
torque exerted on the spheres oscillates quasiperiodically and persistently (af- 
ter an initial transient). The fractional oscillation amplitude (~ 1CU 2 ) increases 
with the angular shear and decreases with the gap width, (ii) When the outer 
\ sphere is accelerated impulsively after an interval of steady differential rotation, 

the torque increases suddenly, relaxes exponentially, then oscillates persistently 
as in (i). The relaxation time-scale is determined principally by the angular ve- 
locity jump, whereas the oscillation amplitude is determined principally by the 



O 

gap width, (iii) When the mutual friction force changes suddenly from Hall- 
Vinen to Gorter-Mellink form, as happens when a rectilinear array of quantized 
Feynman-Onsager vortices is destabilized by a counterflow to form a reconnecting 
vortex tangle, the relaxation time-scale is reduced by a factor of ~ 3 compared 



to (ii), and the system reaches a stationary state where the torque oscillates with 
fractional amplitude ~ 1CU 3 about a constant mean value. Preliminary scalings 
are computed for observable quantities like angular velocity and acceleration as 
functions of Reynolds number, angular shear, and gap width. The results are 
applied to the timing irregularities (e.g., glitches and timing noise) observed in 
radio pulsars. 

Subject headings: dense matter — hydrodynamics — stars: interior — stars: 
neutron — stars: rotation 
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1. Introduction 

The global flow pattern in the superfTuid interior of a neutron star is poorly known, yet 
it is a necessary input into global models of pulsar timing irregularities like glitches (Shemar 
& Lyne 1996; Lyne et al. 2000) and timing noise (Hobbs et al. 2002). The importance of the 
global dynamics was first demonstrated in pioneering laboratory experiments by Tsakadze & 
Tsakadze (1980), who impulsively accelerated rotating spherical and cylindrical containers 
filled with He II and observed abrupt changes in angular velocity, followed by intervals 
of prolonged relaxation, reminiscent of glitches observed in rotation-powered pulsars. The 
relaxation process can be interpreted as Ekman pumping following abrupt spin up (Alpar 
1978; Anderson et al. 1978; Abney & Epstein 1996). The laminar, linear spin up of He II 
between two parallel plates was treated by Reisenegger (1993), who found a relation between 
the poloidal secondary flow and the torque on the container. 

In this paper, we seek to determine which, if any, of the basic features of pulsar tim- 
ing irregularities result directly from the nonlinear hydrodynamics (e.g., oscillations and 
instabilities) of the global flow. To do this, we solve numerically the Hall-Vinen-Bekarevich- 
Khalatnikov (HVBK) equations (Hall & Vinen 1956a,b; Bekarevich & Khalatnikov 1961) 
for a He-II-like superfTuid contained in a differentially rotating spherical shell, building on 
previous numerical simulations of viscous spherical Couette flow in nonastrophysical settings 
(Marcus & Tuckerman 1987a,b; Dumas & Leonard 1994; Mamun & Tuckerman 1995). 

Spherical Couette flow is controlled by two parameters: the aspect ratio (dimensionless 
gap width) S, and the Reynolds number Re. In a slow rotator (Re <C 1), the flow is 
axisymmetric. In a fast rotator (Re ^> 1), the flow is a combination of a primary azimuthal 
rotation and a secondary meridional circulation induced by Ekman pumping (Greenspan 
1968). Transitions from the basic, null vortex state to single and twin vortex states can be 
triggered by gradually increasing Re above a critical value, producing abrupt changes in the 
rotation of the container. This has been observed numerically and experimentally in small 
(5 < 0.24) and large (5 > 0.24) gaps in ordinary fluids (Yavorskaya et al. 1975; Wimmer 
1976; Tuckerman 1983; Sha & Nakabayashi 2001). Hollerbach (1998) also found evidence 
of time-dependent, asymmetric Taylor vortices, in medium gaps (5 = 0.336). However, 
superfluid analogs of these experiments are hard to perform, partly because of problems 
with visualization techniques (Bielert & Stamm 1993). Numerical HVBK simulations have 
been performed in cylindrical Taylor- Couette geometries (Henderson et al. 1995; Henderson 
2001; Henderson & Barenghi 2000), delivering good agreement with experiments, but not 
yet in the more challenging spherical Couette geometry. We attempt this here for the first 
time. 

The paper is organized as follows. In §2, we establish the astrophysical context of our 
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numerical experiments by explaining how the idealized spherical Couette problem relates to 
a realistic neutron star. HVBK theory is reviewed in §3. The numerical method in spherical 
geometry is described and verified against test cases in §4. Results are then presented for the 
long-term response of the superfluid in three scenarios: the outer and inner spheres rotate 
steadily and differentially (§5), the outer sphere is accelerated impulsively (§6), and the 
mutual friction force changes suddenly (from anisotropic to isotropic form) in response to a 
line vortex instability (§7). Finally, in §8, we apply the results to pulsar timing irregularities 
and discuss the limitations of our approach. 



2. Idealised hydrodynamic model of the outer core of a neutron star 

In this paper, we conduct numerical experiments that seek to mimic, in an idealised 
context, the flow of neutron superfluid within the outer core of a rotating neutron star. The 
density of the outer core lies in the range 0.6p* < p < 1.5p*, where p* = 2.6 x 10 14 g cm~ 3 
is the nuclear saturation density, which translates to radii in the range 5 km <, r <, 9 km 
(Sedrakian & Sedrakian 1995; Weber 1999; Yakovlev et al. 1999; Dean & Hjorth- Jensen 
2003). We focus on the outer core because the superfluid is probably unpinned in this 
region (see below) — that is, it behaves essentially hydrodynamically — and its physical 
state is reasonably well understood, unlike the inner core. In the upper part of the outer 
core, 0.6p* < p < p*, the neutrons form Cooper pairs in the 1 5'o state; in the lower part, 
p* < p < 1.5p*, they pair in the 3 P2 state (Weber 1999). Here, for simplicity, we assume that 
the isotropic 1 S'o phase fills the entire outer core. The 3 P2 phase is anisotropic: gradients in 
the orientation (texture) of the pair angular momenta drive counterflows between the viscous 
and inviscid components of the superfluid (Vollhard & Wolfle 2002; Mastrano & Melatos 
2005), introducing complications (e.g. boundary conditions at the 1 S - 3 P 2 interface) which 
lie outside the scope of this paper. 4 

The outer core contains charged species, particularly protons and electrons, which are 
incorporated into the viscous component of the superfluid in our simulations. The protons 
are probably in a type II superconducting state, but the electrons are not (Sauls 1989). 
In a type II superconductor, the protons inside the magnetic fluxoids interact with the 
Feynman-Onsager vortices in the neutron superfluid, serving as natural pinning sites (Sauls 
1989); the pinning geometry can be complicated (Ruderman et al. 1998) if the core magnetic 



4 A two-stream (Andersson et al. 2003) or Kelvin-Helmholtz (Mastrano & Melatos 2005) instability can 
be excited at the interface if the isotropic and anisotropic phases are in relative motion. 
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field has comparable poloidal and toroidal components (Thompson & Duncan 1993). 5 
We neglect this effect in what follows, as well as the effects arising from entropy gradients 
(Mastrano & Melatos 2005). In addition, the protons exert an effective drag on the neutrons 
by the so-called entrainment effect (Mendell 1991; Sedrakian & Sedrakian 1995), in which 
the momentum of one species is partly carried along by the other species, as in 3 He- 4 He 
mixtures (Andreev & Bashkin 1976; Andersson & Comer 2001). We neglect this effect in 
the computations in this paper, where there is no entropy gradient at the 3 P 2 - 1 S , interface 
(Mastrano & Melatos 2005). 

The hydrodynamic boundary conditions at the upper and lower surfaces of the outer 
core are set by the conditions in the inner crust and inner core respectively. The physical 
state of the inner core is essentially unknown, so we experiment with the two extremes of 
rough walls (pinning) and perfect slip (no pinning) at this boundary. For example, pinning 
is favored if triple-flavor color superconductivity results in a crystal lattice in the inner core 
(Rajagopal 2002). In the inner crust, it is widely believed that the Feynman-Onsager vortices 
are trapped at nuclear pinning sites (Sauls 1989), so that a rough wall is the appropriate 
boundary condition. However, recent work suggests that, even in the inner crust, pinning 
may not occur. Donati & Pizzochero (2003) studied the vortex-nucleus interaction and found 
intersticial pinning only, with low pinning forces ^ 0.4 MeV fm _1 . 

In the standard picture of pulsar glitches, the stellar crust is loosely coupled to a rapidly 
rotating superfluid interior (Lyne & Graham-Smith 1998). Although we do not attempt to 
explicitly construct a model of glitches in this paper, our simulations address a similar 
scenario: the upper and lower surfaces of the outer core rotate differentially, as the crust 
spins down electromagnetically, until a glitch occurs and the crust suddenly accelerates. We 
study the global superfluid flow pattern in the outer core before and after the glitch. We also 
study how the Feynman-Onsager vortices switch from a rectilinear array to a reconnecting 
tangle via the Donnelly-Glaberson instability (Glaberson et al. 1974). The transition from 
an organized to a disorganized distribution of superfluid vorticity in the outer core, together 
with the acceleration of the crust, strongly affects the rotation of the star, as we show below. 



5 Pinning in the outer core is arguably inconsistent with the precession periods (~ 1 yr) observed in some 
pulsars (Stairs et al. 2000), because the Feynman-Onsager vortices damp the precession on short time scales 
(~ 1 hr) as they pass through the array of magnetic fluxoids (Link 2003). 
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3. HVBK theory 

HVBK theory is a generalization of the two-fluid Landau-Tisza theory for superfluid He 
II that includes the physics of quantized Feynman-Onsager vortices (Hall & Vinen 1956a,b; 
Bekarevich & Khalatnikov 1961). It applies equally well to a neutron superfluid with 1 J S , 
parity (Tilley & Tilley 1986). Fluid particles in the theory are assumed to be threaded by 
many coaligned vortices. In the continuum limit, we can define a macroscopic superfluid 
vorticity w s = Vxv s which is not zero, despite the fact that, microscopically, the superfluid 
obeys V x v s = . This is valid if the length scales in the flow are longer than the average 
separation between vortex lines. The isothermal HVBK equations of motion take the form 
(Barenghi & Jones 1988; Henderson & Barenghi 2000) 

— — = + z/„Vv„ + — F Vw,, (1) 

at p p p 



^^^^.T.ftlF-^Vkl, (2) 
at p p p 

with d n , s /dt = d/dt + v„ iS • V, supplemented by the incompressibility condition V • v n = 
V- v s = 0, a good approximation in a neutron star, where the flow is subsonic. In (1) and (2), 
v njS and p ntS are the normal fluid and superfluid velocities and densities respectively. Effective 
pressures p s and p n are defined by Vp s = Vp — |/0nV(v^ s ) and Vp n = Vp + §p s V(v£ s ), 
with v ns = v„ — v s . We define v n to be the kinematic viscosity of the normal fluid and 
v s = (k/Att) log(bo/a ) to be the stiffness parameter, where k = h/2m n is the quantum of 
circulation, m n is the mass of the neutron, ao is the radius of the vortex core, and bo is 
the intervortex spacing. Note that v s has the same dimensions as a kinematic viscosity but 
controls the oscillation frequency of Kelvin waves excited on the vortex lines (Henderson 
et al. 1995). 

The vortex tension force per unit mass, v s T, is given by 

T = w s x(Vx u 8 ), (3) 

with lj s = u; s /|u; s |. The tension arises from local circulation around quantized vortex lines. 
It corresponds to the Magnus self-force exerted on a line by its own, self-induced velocity 
field, which is given by u s u s x (V x Cj s ) in the local induction approximation of the full 
Biot-Savart integral, when the radius of curvature is much larger than a (Donnelly 1991). 

The mutual friction force per unit mass arises from the interaction between the quantized 
vortex lines and the normal fluid (via roton scattering in He II and electron scattering in 
a neutron star) (Hall & Vinen 1956a,b). The form of this force depends on the global 
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configuration of the vortices. If the vortices ccupy a rectilinear array, the friction force per 
unit mass takes the Hall-Vinen (HV) form (Hall & Vinen 1956a,b) 

F = ^Bu s x (u a x v ns - u s T) + )rB'(u s x v ns - u s T) (4) 

where B and B' are dimensionless, temperature-dependent coefficients (Barenghi et al. 1983). 
However, in a realistic neutron star, the distribution of vorticity need not correspond to a 
rectilinear vortex array (Greenstein 1970, 1974). For example, differential rotation can trigger 
the Donnelly-Glaberson (DG) instability (Cheng et al. 1973; Glaberson et al. 1974), in which 
a dense tangle of vortices forms in the counterflow. Experimental evidence for this effect 
comes from counterflow experiments that measure the attenuation of second sound in narrow 
channels (Vinen 1957; Swanson et al. 1983), supported by numerical simulations using the 
vortex filament method (Tsubota et al. 2003). From these experiments, it is inferred that 
the friction force per unit mass takes the isotropic Gorter-Mellink (GM) form (Gorter & 
Mellink 1949) 

F = A , (Wlj Vns , (5) 

where A' = B 3 p 2 l n 2 Xi/3p 2 X2 is a dimensionless, temperature-dependent coefficient, related 
to the original GM constant (usually denoted as A in the literature) by A' = Apn. Here, 
Xi and X2 are dimensionless constants of order unity (Vinen 1957). Note that (1) and (2), 
derived assuming a dense, roughly parallel array of vortices, may not be consistent with (5). 
The validity of HVBK theory in turbulent flow is yet to be established (Barenghi et al. 1995). 



4. Numerical Method 

4.1. Pseudospectral collocation in spherical geometry 

Equations (1) and (2) are discretized in space using a pseudospectral collocation method 
(Canuto et al. 1988). A time-split fractional step algorithm advances the solution (Canuto 
et al. 1988; Bell & Colella 1989; Streett & Hussaini 1991). Nonlinear terms are treated explic- 
itly using a third-order Adams-Bashforth scheme and diffusive terms are treated implicitly 
using a Crank-Nicholson scheme (Boyd 2001). Incompressibility is enforced by pressure- 
correction projection (Brown et al. 2001; Chorin 1968). We follow closely the approach of 
Bagchi & Balachandar (2002), expanding the polar (9) and azimuthal (0) coordinates in 
Fourier functions and the radial (r) coordinate in Chebyshev polynomials. The velocity 
fields are expressed as u(r,9,(f)) = J2i j k ( ~'ijkTi(r)f(j9)e' lk 't > , where Tj(r) is the i-th Cheby- 
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shev polynomial, f(8) equals sin# or cos6*, and i, j, k are integers running from to N r , Ng, 
respectively. The pole parity problem which plagues Fourier series on a sphere (Boyd 
2001) must be handled with care; the tangential and azimuthal components of v n and v s 
change sign across the poles (Orszag 1974; Fornberg 1998). 

Spectral methods are global and therefore sensitive to boundary conditions. The normal 
fluid v n satisfies no-slip and no-penetration boundary conditions at the inner and outer 
surfaces of the container. There is no general agreement in the literature on what boundary 
conditions are suitable for v s . Khalatnikov (1965) suggested that vortex lines can either 
slide along, or pin to, the boundaries. The tangential velocity v^, of a vortex line relative to 
a rough boundary is given by (Hills & Roberts 1983) 



where n is the unit normal to the surface; in this case, vortex lines are permanently attached 
to the surface. At a smooth boundary, on the other hand, one has 



in this case, vortex lines are oriented perpendicular to the surface. Condition (6) is difficult 
to implement in HVBK theory, in which each fluid element is threaded by many vortex lines, 
because v L cannot be calculated from v n and v s . On the other hand, frictionless sliding 
(u s x n = 0) leads to numerical instabilities, because u s can develop parallel components 
near the surface as it evolves (Henderson & Barenghi 2000). Therefore, in this paper, we 
adopt no-slip and no-penetration [i.e., pinning; see Baym et al. (1992)] boundary conditions 
for v s to stabilize the evolution. 

Spectral methods can develop global oscillations, due to the Gibbs phenomenon, when 
applied to problems with stiff numerical solutions or discontinuities (Vandeven 1991). In the 
superfluid, such oscillations are damped indirectly (and weakly) through F. To mitigate this 
problem, we filter out high spatial frequency modes in the r and 9 expansions by multiplying 
the spectral coefficients by an exponential filter defined by a e {k/N r ^) = exp[— (k/N r ^ e ) y lne], 
with < \k\ < N r fi, where e = 2.2 x 10~ 16 is the machine zero and 7 is the order of the 
filter. Strong filtering is needed (7 < 8) in order to stabilize the evolution. Most of our runs 
employ grids with (jV r , N , iV^) = (210, 250, 4), although a few have = 32. 

As our simulations are the first to treat the spherical Couette problem for superfluids, 
we are obliged to validate the code in the simpler case of a classical viscous fluid. In this limit 
(p s — > 0, p n — > p as T — > T c ), we succesfully reproduce the — >• 1 and 1 — > 2 vortex-state 
transitions observed in small and large gaps in spherical Couette flow (Marcus & Tuckerman 



(n • u) s )u s x v L = 0, 



(6) 



(n ■ uj s )uj s x n = 0; 



(7) 
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1987a,b; Dumas & Leonard 1994), achieving agreement to three significant digits in the 
critical value of Re and nine significant digits in the torque. 

4.2. Model pulsar 

Our model pulsar consists of a spherical shell filled with a 1 S , -paired neutron superfluid. 
The inner (radius Ri) and outer (radius R 2 ) surfaces of the shell rotate at angular velocities 
Qi and Q 2 , respectively, about a common axis. Henceforth, all quantities are expressed in 
dimensionless form using R 2 as the unit of length and Q^ 1 as the unit of time. The Reynolds 
number and dimensionless gap width are defined by Re = ^l\R\jv n and 5 = (R 2 — Ri)/R 2 
respectively. 

Our choice of parameters is limited by computational capacity. For example, in a 
realistic pulsar, one might expect 10~ 9 ^ AQ/Q ^ 10~ 6 in the lead-up to a glitch (Hobbs 
2002), with AQ = \Q 2 — We take 0.1 < AQ < 0.3 instead, in order to observe the 
long-term effect of differential rotation (i.e., several full "wraps" of the inner sphere relative 
to the outer sphere) over a time interval < t < 90 that is computationally achievable. 

Initially, both fluids must satisfy the incompressibility condition, with (v ra ) r = (v s ) r = 
at r = Ri, R 2 . We choose the Stokes solution (Landau & Lifshitz 1963) as the initial 
condition, with v s = v n . We also take p s /p = 0.908 and p n /p = 0.092, corresponding to 
the superfluid and normal fluid fractions in He II at a fiducial temperature T/T c = 0.66 
(in a realistic pulsar, one has T/T c ^ 0.1). Importantly, the Reynolds number in a typical 
neutron star can reach Re ~ 10 11 (Mastrano & Melatos 2005), whereas, in our simulations, 
we restrict ourselves to 10 2 < Re < 10 4 for two reasons: for Re > 10 5 , the normal fluid flow 
becomes turbulent, and we cannot resolve it properly; and, when Re is smaller, a steady 
state is reached in a shorter time. 

We set the stiffness parameter to v s = 10~ 5 cm 2 s _1 , such that v s <^ v n as in a young 
pulsar (Greenstein 1970), even though a more realistic value is v s = 10" 3 cm 2 s _1 [cf. v n 
10 2 cm 2 s _1 ; Mastrano & Melatos (2005)]. In practice, v s cannot be chosen independently of 
5. For a small gap (5 < 0.1) or large z/ s , the tension force disrupts the flow and leads to 
numerical instabilities, in which the stiff superfluid streamlines bend sharply at the walls. 
This effect is suppressed in a cylindrical shell, where the tension causes the superfluid to 
rotate as a column parallel to the curved walls (Henderson & Barenghi 2000). Consequently, 
perfect slip boundary conditions on v s may prove beneficial in the spherical problem; we will 
investigate them in a future paper. 

The HV friction force parameters B and B' are unknown in a neutron star. We choose 



- 9- 



B = 1.35 and B' = 0.38, adopting the He II values at T/T c = 0.66 (Barenghi et al. 1983; 
Donnelly 1991; Donnelly & Barenghi 1998). The parameter A' = 5.8 x 10~ 3 (with X1/X2 — 
0.3) at the same temperature can be calculated from a fitting formula derived by Dimotakis 
(1972), which is consistent with previously published experimental values (Vinen 1957). 
Stable long-term evolution is difficult to achieve for this value of A', so we take A' = 5.8 x 10~ 2 
instead. A more detailed study of the dynamics for a wider range of B, B' and A' will be 
presented in a forthcoming paper. 

5. Steady differential rotation 

In this section, we describe the results of a control experiment in which the the inner 
and outer spheres rotate steadily yet differentially, i.e., Qi and Q2 7^ ^1 are held constant. 
We consider four cases, identified by Au, Bu, Cu, and Du in Table 1, corresponding to 
medium and large gaps (5 = 0.5, 0.3) and small and medium shear (0.7 < f2 2 < 0.9). A few 
preliminary runs with N$ = 32 confirm that the flow remains axisymmetric, so we restrict 
the resolution to (N r , No, N^) = (150,400,4) in what follows. 

Figure 1 displays the meridional streamlines of the normal fluid and the superfluid for 
5 = 0.3 and AQ = 0.3 as a function of time. The corresponding torques on the inner 
and outer spheres are plotted as dashed curves in Figures 2a and 2b respectively, with 
squares marking the instants when the streamlines are plotted in Figure 1. The torque 
oscillates quasi-periodically, with a final periodic state that persists as long as the differential 
rotation is maintained (up to t = 90 in our runs), with a peak-to-peak fractional amplitude 
AN Z /N Z m 0.02. The flow pattern in the normal fluid contains a secondary vortex near 
the poles for 8 < t < 20, which expands and contracts quasi-periodically. For t > 20, 
the secondary vortex disappears and a single meridional shell remains, which also expands 
and contracts, but only slightly. The superfluid streamlines display a richer pattern: an 
equatorial vortex develops at 6 < t < 8, a polar vortex develops at 8 < t < 20, and the 
mid-latitude flow becomes disorganized for t > 20. Note that we must reduce the spatial 
resolution to (N r , N$) = (120,250) during the run in order to stabilize the evolution; for 
(N r ,Ng) = (150,400), one observes sudden jumps in the torque unless At < 10~ 5 (which 
lengthens the runs unacceptably). It is possible that turbulence is generated in v s for t > 20 
and that, by intentionally choosing a coarser resolution, we are effectively smoothing over 
the turbulent eddies. 

The torque evolves similarly in all four cases. Initially, the linear Stokes solution adjusts 
to a nonlinear, cellular solution (with meridional circulation) within a time At ~ 6. This 
transient state is of no interest astrophysically. Subsequently, the torque oscillates persis- 
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tently — a new and astrophysically relevant phenomenon. The oscillation amplitude depends 
strongly on the gap width 5 and weakly on the shear Afl. For example, if 5 = 0.5 is held 
constant, the amplitude decreases from AN Z /N Z = 0.026 to AN Z /N Z = 0.014 as the shear 
decreases from AQ = 0.3 (solid curve in Figure 2) to AQ = 0.1 (dotted curve in Figure 2). 
On the other hand, if AQ = 0.3 is held constant, the amplitude decreases from 0.036 to 
0.026 as S increases from 0.3 to 0.5 (dashed and solid curves in Figure 2). 

6. Impulsive acceleration of the outer sphere 

When a glitch is triggered in a rotation-powered pulsar, the crust and consequently the 
superfluid interior are spin up impulsively. We simulate a spin-up event of this kind by 
abruptly accelerating the rotation of the outer sphere from 0.7 < ^2 < 0.9 to Q2 = 1 at 
t — 20. Figure 3a shows the change in the flow pattern of the normal fluid and superfluid 
before and after this sudden acceleration. The difference is more pronounced in the normal 
fluid; a secondary vortex appears at mid-latitudes in each hemisphere, adjacent to the outer 
shell. The jump in torque comes almost entirely from these mid-latitude regions. 

Figure lb shows how the dominant spectral coefficients (k = 1) of (v n )# behave 
before and after the spin- up event. The three modes of highest amplitude before acceler- 
ation are C221 > C 42 i > C 52 i- After acceleration, the same ordering is observed, but with 
amplitudes 12% lower, and the fourth-strongest mode switches from C 62 i to C 32 i. By 
contrast, in similar data for (v s )$, no change is seen in the ordering of the top 10 spectral co- 
efficients, and their amplitudes change by less than 10%. The superfluid is therefore largely 
unaffected by the glitch, as one might anticipate by inspecting Figure 3; most of the effect 
is transmitted to the normal fluid by the viscous Ekman layer adjacent to the outer shell. 

A key numerical issue is whether the spectral method faithfully resolves the flow. Em- 
pirically, this occurs if the mode amplitudes decrease quasi-monotonically with polynomial 
index. We confirm this here. The Chebyshev modes, plotted in Figure 4, decrease by a 
factor ~ 10 4 in amplitude from % — 1 to % — 50, indicating that the flow is resolved in r. 
We also increase the resolution in at t > 20 by Fourier interpolation and confirm that the 
results are identical for 4 < N$ < 32. 

When the outer sphere is accelerated, we find that the time-step must be halved to avoid 
sudden artificial jumps in the torque. This property has been noted widely in the literature: 
in viscous spherical Taylor- Couette flow, the nonlinear dynamics depends sensitively on the 
time-step in the transition between different vortex states (Loukopoulos & Karahalios 2004; 
Sha & Nakabayashi 2001; Liu et al. 1996). For example, Loukopoulos & Karahalios (2004) 
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reported that the transition to an asymmetric 1-vortex state can only be achieved by varying 
Re quasi-statically near its critical value and using a small (At = 1CT 4 ) time step. Different 
flow states in viscous spherical Couette flow are associated with energy differences. If the 
system is driven through an intermediate state too quickly, the flow cannot follow the fast 
change, due to its inertia, and a steady final state cannot be reached (Wimmer 1976). 

The post-acceleration evolution of the torque is qualitatively similar in the four cases 
studied. This can be observed in Figure 5, where we plot the fractional change in the angular 
velocity AQ/Q and its first time derivative Ati/Cl as a function of time. 6 In radio pulsars, 
both these quantities are observed to high accuracy in radio timing experiments. The outer 
torque jumps at t — 20, oscillates with a period At ~ 5, then decays quasi-exponentially. 
Note that, in Figure 56, the curves for Aa, Ca and Da are grouped together quite closely 
before acceleration, whereas the curves for Aa and Ba are grouped closely after acceleration. 
This implies that the gap width 5 is the main parameter controling the spin-down torque 
over the long term (i.e., during steady differential rotation), as in §5, whereas the angular 
velocity jump |f^2 — ^i| is the main parameter controling the relaxation time-scale. By 
fitting an exponential to the secular evolution of AQ, we find that the (e -1 ) relaxation time- 
scale is given by r = 12, 9.5, 9.3, and 8.2 for Aa, Ba, Ca, and Da respectively and scales 
roughly as r oc <5 1//3 1 — ^i| 1//3 , similar to (but not exactly the same as) the Ekman scaling 
roc l^-fiil 1 / 2 . 



7. Instability of the Vortex Array 

The stability of a rectilinear array of quantized vortices was investigated by Vinen in 
counterflow experiments in narrow channels (Vinen 1957). He attributed the attenuation of 
second sound to the generation of a "vortex tangle" in the counterflow. The critical (axial) 
counterflow velocity v ns for the onset of this turbulent state was calculated theoretically 
by Glaberson et al. (1974), for small perturbations of a single vortex line (Cheng et al. 
1973). He found that, for v ns > 2(2f2z/ s ) 1 / 2 , growing Kelvin waves are excited and the 
vortex lines reconnect to form a tangle, in accord with experimental data (Cheng et al. 
1973). Similar results have been produced by numerical simulations with the vortex filament 
method (Schwarz 1988; Tsubota et al. 2003). The inclusion of slow rotation in counterflow 
experiments reduces the critical v ns at which superfluid turbulence sets in (Swanson et al. 



6 Note that Sl2(t) is held fixed in the code after the acceleration at t = 20. Hence the angular velocity 
change Af2 = 2 (i) — O 2 (20) plotted in Figure 5 is not computed self-consistently. It is found by integrating 
the viscous torque on the outer sphere, calculated by solving (1) and (2) subject to the time-independent 
boundary condition ri 2 (t) = constant. 
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1983). 

In order to simulate what happens when a vortex tangle develops from an initially 
regular vorticity distribution, we suddenly change the friction force in (1) and (2) from HV 
to GM form, according to the prescription in §4. From Figure 6, we observe that the torque 
reaches a stationary state ~ 3 times more quickly than in the spin-up experiments presented 
in §6. The torque decays more quickly than an exponential, as can be appreciated from 
Figure 6. Nevertheless, we can estimate the (e _1 ) relaxation time-scale for the cases Af, 
Bf, Cf, and Df in Table 1 to be 3.4, 3.4, 3.0, and 3.6, respectively. The relaxation time is 
essentially independent of 5 and |f2 2 — 

Relaxation occurs more quickly than in the spin-up experiment because the GM force is 
small compared to the HV force. From §6, we have |F GM |/|F H v| = (A' '/ ' B)(p n p s / p 2 ){v^ s / 'kuj s ) 
[and hence oc (p n / p) 3 for p n ^ p\; on the other hand, we find |Fqm|/|Fhv| ~ 10 -5 empirically 
from our runs, implying v 2 s /ku s ~ 2.2 x 10~ 6 . The sudden change in the mutual friction 
effectively uncouples the normal fluid and the superfluid, and the torque evolves rapidly to 
a constant value on a time-scale At ~ 3. After this time, the inner and outer torques are 
equal to one part in ~ 10 3 and their mean values hardly vary at all. Nevertheless, while 
the torque appears constant for t > 25 in Figure 6, closer examination on a magnified scale 
shows that \N 2 — N±\ oscillates persistently, out to t = 90, with peak-to-peak amplitude 
AN Z /N Z ~ 10~ 3 . 



8. Discussion 

In this paper, we present the first three-dimensional hydrodynamic HVBK simulations 
of a 1 So-paired neutron superfluid in a rotating spherical shell, generalizing studies of viscous 
fluids in spherical Couette geometry (Marcus & Tuckerman 1987a,b; Dumas & Leonard 1994) 
and He II in cylindrical Couette geometry (Henderson et al. 1995; Henderson & Barenghi 
2000). The code accurately resolves the superfluid flow pattern for Reynolds numbers and 
gap widths in the range 10 < Re < 10 4 and 0.2 < 5 < 0.7, respectively. 

Persistent quasiperiodic oscillations are always observed in the torque during steady 
differential rotation, after the initial transient dies away, with typical period ~ fi^ 1 and 
fractional amplitude ~ 10~ 2 . The oscillation amplitude increases as Re increases. The 
gap width 5 principally controls the spin-down torque over the long term, during steady 
differential rotation, whereas the frequency jump \Q 2 — &i\ principally sets the relaxation 
time-scale after a sudden acceleration of the outer sphere. If, instead of spinning up the 
outer sphere, we suddenly switch the mutual friction force from anistropic (HV) form to 
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isotropic (GM) form, as happens when a rectilinear array of quantized vortices reconnects 
unstably to form a tangle, the relaxation time-scale is ~ 3 times shorter and the torque 
difference \N 2 — N ± \ oscillates subsequently with peak-to-peak fractional amplitude ~ 1CT 3 
around a constant mean value. We speculate that a transition from a turbulent flow, driven 
by long-term differential rotation, to a laminar flow immediately after a glitch, might cause 
some of the timing irregularities observed in pulsars. Similar transitions from turbulent to 
laminar flow in a superfluid have already been observed in laboratory experiments where 
He II, cooled to a few mK, flows around an oscillating microsphere (Niemetz et al. 2002; 
Schoepe 2004). 

Computational limitations impose upon us some artificial approximations. Some of 
these approximations will be tested independently in a future paper. For example, one 
might hope to use a weaker filter when damping the high spatial frequency modes, perhaps 
by introducing an artificial viscosity (Andreassen et al. 1994; Ould Kaber 1996), and one 
can model small-scale turbulence using the techniques of large eddy simulations (Germano 
et al. 1991; Lilly 1992; Mahesh et al. 2004). Moreover, although some of our parameters are 
not as close as one might which to realistic neutron star values, we are careful to respect 
the ordering of physical quantities encountered in a neutron star (e.g., Af2 2 ^ ^2) Re ^ 1, 
v s <C v n ) and to use values of the dimensionless friction coefficients (B, B' and A') that 
apply to the neutron star regime p n <C p s and are consistent with experimental data on He 
II. By artificially reducing u s by two orders of magnitude, we trade off some accuracy for 
computational speed in our modeling of old pulsars, where the tension force is the primary 
interaction (Greenstein 1970), although our modeling of young pulsars is mostly unimpaired 
(Hobbs et al. 2004; Lyne et al. 2000). 

We consider superfluid flow in the outer core of the star, where there is no pinning 
in the body of the fluid, only at the boundary. Consequently, those parts of our analysis 
pertaining to the formation of a reconnecting vortex tangle via the DG instability may not 
translate straightforwardly to the inner crust, where pinning at multiple lattice sites along 
a vortex line may suppress the DG instability and even prevent a tangle from forming at 
all. On the other hand, if Ekman pumping is the primary core-crust coupling mechanism in 
a real neutron star, then the dynamics in the outer core must dominate the physics; crust- 
only mechanisms are inconsistent with observations of Vela's glitches in the Ekman regime 
(McCulloch et al. 1990; Abney et al. 1996). 

The time-scales of torque fluctuations in our numerical experiments are generally ^ lOf^ 
much shorter than the time-scales for post-glitch relaxation and timing noise oscillations ob- 
served in radio pulsars. However, it is crucial to realize that the numerical time-scales will 
lengthen greatly in more realistic simulations, where the angular shear is lower and Re is 
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higher. Such numerical experiments are computationally challenging, and we are investigat- 
ing alternative numerical approaches that may render them tractable. 

A necessary extension of the code is to include radial stratification. It has been argued 
that, for a viscous fluid, stratification hampers the formation of an Ekman layer near the 
outer boundary (Abney & Epstein 1996). The Brunt- Vaisala frequency in a neutron star was 
estimated to be ubv ~ 500s _1 around nuclear density and cjbv ~ 0.7^2 for the Vela pulsar 
(Reisenegger & Goldreich 1992), whereas one needs uj-by ~ 0.2Q 2 for Ekman pumping to pro- 
ceed effectively (Abney & Epstein 1996). Stratification can also affect the coupling between 
the crust and the superfluid core of a neutron star by modifying precessional modes (Levin 
& D'Angelo 2004). The pressure projection operation in our code relies on incompressibility, 
so stratification is hard to include self-consistently. As a first approximation, however, one 
can set v r = in the incompressible calculation, as proposed by Levin & D'Angelo (2004). 

The HVBK equations may not apply to a vortex tangle. An alternative theory which 
explicitly treats displacements and oscillations of the vortex array was developed by Chandler 
& Baym (1983, 1986). The complexities thereby introduced, and fresh uncertaintities over 
boundary conditions, raise new numerical challenges. 

We acknowledge the computer time supplied by the Australian Partnership for Ad- 
vanced Computation (APAC) and the Victorian Partnership for Advanced Computation 
(VPAC). We also thank George Hobbs for helpful discussions and Caroline Andrzejewski 
for undertaking a thorough literature review of spherical Couette flow in support of this 
research. 
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Table 1: Simulation Parameters 
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Fig. 1. — Meridional streamlines for the normal fluid (left) and the superfluid (right) as a 
function of time (6 < t < 50) for 5 = 0.3 and AQ = 0.3 (case Bu, in Table 1). The time is 
expressed in units of . 
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Fig. 2. — (a) Evolution of the (a) inner torque and (b) outer torque for the cases Au (solid 
curve), Bu (dashed curve), Cu (dash-dotted curve), and Du (dotted curve), whose parameters 
are quoted in Table 1. The time is expressed in units of il^ 1 and the torques in units of 
pR\Sl\. The filled squares correspond to the nine instantaneous snapshots of the streamlines 
plotted in Figure 1. 




Fig. 3. — (a) Meridional streamlines before (left pair of hemispheres) and after (right pair 
of hemispheres) the outer sphere is accelerated instantaneously from Q 2 — 0.7 to f2 2 = 1 at 
t = 20. For each pair of hemispheres, the normal fluid streamlines are plotted at left and the 
superfluid streamlines are plotted at right. The results are for the case Af in Table 1. (b) 
Histograms of the magnitudes of the 12 dominant spectral coefficients of (v n )o before (left) 
and after (right) the acceleration at t = 20, as a function of the Chebyshev polynomial index 
i and the polar Fourier index j, with azimuthal index k — 1. 
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Fig. 4. — Modes amplitudes for the r (o), 9 (□), and (*) components of the normal fluid 
and the r (|), 6 (x), and <fi ((}) components of the superfluid as a function of the Chebyshev 
polynomial index i, at t — 20.1. The results are for case Af in Table 1. 
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Fig. 5. — (a) Fractional change in angular acceleration of the outer sphere, AQ/Q = [f2 2 (^) ~~ 
fi 2 (20)]/fi 2 (20), as a function of time, (b) Fractional change in angular velocity of the outer 
sphere, Ail/il = [f2 2 (*) — ^2(20)]/fi 2 (20), as a function of time. Time is measured in units 
of n^" 1 . Four cases are shown, whose parameters are quoted in Table 1: Aa (solid curve), Ba 
(dashed curve), Ca (dashed-dotted curve), and Da (dotted curve). 




time time 

Fig. 6. — (a) Fractional change in angular acceleration of the outer sphere, Ati/Cl = [il 2 (£) — 
^2 (20)] /^2 (20), as a function of time, (b) Fractional change in angular velocity of the outer 
sphere, Ail/il = [Q 2 (t) — il 2 (20)]/fi 2 (20), as a function of time. Time is measured in units 
of Jlj" 1 . Four cases are shown, whose parameters are quoted in Table 1: Af (solid curve), Bf 
(dashed curve), Cf (dashed-dotted curve), and Df (dotted curve). 



